# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating figures and tables in Appendix D
# Appendix D Table 8: Influence of Logged Crime on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix D Table 9: Influence of Logged Crime on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status.
# Appendix D Figure 6: Louisville: Influence of Logged Crime on Logged Arrests By Type (Standardized), Conditional on Racial Boundary Status.
# Appendix D Table 10: Louisville: Influence of Logged Crime on Logged Arrests by Type (Standardized), Conditional on Racial Boundary Status.



suppressPackageStartupMessages(
  
  {
    library(AER)
    library(dplyr)
    library(fixest)
    library(lfe)
    library(tidyverse)
    library(ggplot2)
    library(MASS)
    library(sensemakr)
    library(haven)
    library(readstata13)
    library(readxl)
    library(readr)
    library(gridExtra)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggthemes)
    library(meta)
  }
)

# first prepare variables by city
# read in data by city

# Load Atlanta final dataset
load("atl_final.RData")

# log and +1 to variables (pop already logged)
atl_fin = atl_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
atl_fin$larrest_sd <- atl_fin$larrests - (mean(atl_fin$larrests)/sd(atl_fin$larrests))
atl_fin$lmisdemeanors_sd <- atl_fin$lmisdemeanors - (mean(atl_fin$lmisdemeanors)/sd(atl_fin$lmisdemeanors))
atl_fin$lfelonies_sd <- atl_fin$lfelonies - (mean(atl_fin$lfelonies)/sd(atl_fin$lfelonies))
atl_fin$lnonviolent_sd <- atl_fin$lnonviolent - (mean(atl_fin$lnonviolent)/sd(atl_fin$lnonviolent))
atl_fin$lviolent_sd <- atl_fin$lviolent - (mean(atl_fin$lviolent)/sd(atl_fin$lviolent))
atl_fin$lsociety_sd <- atl_fin$lsociety - (mean(atl_fin$lsociety)/sd(atl_fin$lsociety))
atl_fin$lperson_sd <- atl_fin$lperson - (mean(atl_fin$lperson)/sd(atl_fin$lperson))
atl_fin$lproperty_sd <- atl_fin$lproperty - (mean(atl_fin$lproperty)/sd(atl_fin$lproperty))

# rename racial and economic boundary variable
atl_fin$atl_white_blv <- atl_fin$p_race_white_blv
atl_fin$atl_ses_blv <- atl_fin$ses_blv

# Load Austin final dataset
load("aus_final.RData")



# log and +1 to variables (pop already logged)
aus_fin = aus_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_fin$larrest_sd <- aus_fin$larrests - (mean(aus_fin$larrests)/sd(aus_fin$larrests))
aus_fin$lmisdemeanors_sd <- aus_fin$lmisdemeanors - (mean(aus_fin$lmisdemeanors)/sd(aus_fin$lmisdemeanors))
aus_fin$lfelonies_sd <- aus_fin$lfelonies - (mean(aus_fin$lfelonies)/sd(aus_fin$lfelonies))
aus_fin$lnonviolent_sd <- aus_fin$lnonviolent - (mean(aus_fin$lnonviolent)/sd(aus_fin$lnonviolent))
aus_fin$lviolent_sd <- aus_fin$lviolent - (mean(aus_fin$lviolent)/sd(aus_fin$lviolent))
aus_fin$lsociety_sd <- aus_fin$lsociety - (mean(aus_fin$lsociety)/sd(aus_fin$lsociety))
aus_fin$lperson_sd <- aus_fin$lperson - (mean(aus_fin$lperson)/sd(aus_fin$lperson))
aus_fin$lproperty_sd <- aus_fin$lproperty - (mean(aus_fin$lproperty)/sd(aus_fin$lproperty))

# rename racial and economic boundary variable
aus_fin$aus_white_blv <- aus_fin$p_race_white_blv
aus_fin$aus_ses_blv <- aus_fin$ses_blv


## Load Boston data
load("bos_final.RData")


# log and +1 to variables (pop already logged)
bos_fin = bos_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime + 1))



# create scaled DVs (to account for different pop. sizes and density in blocks)
bos_fin$larrests_sd <- bos_fin$larrests - (mean(bos_fin$larrests)/sd(bos_fin$larrests))
bos_fin$lmisdemeanors_sd <- bos_fin$lmisdemeanors - (mean(bos_fin$lmisdemeanors)/sd(bos_fin$lmisdemeanors))
bos_fin$lfelonies_sd <- bos_fin$lfelonies - (mean(bos_fin$lfelonies)/sd(bos_fin$lfelonies))
bos_fin$lnonviolent_sd <- bos_fin$lnonviolent - (mean(bos_fin$lnonviolent)/sd(bos_fin$lnonviolent))
bos_fin$lviolent_sd <- bos_fin$lviolent - (mean(bos_fin$lviolent)/sd(bos_fin$lviolent))
bos_fin$lsociety_sd <- bos_fin$lsociety - (mean(bos_fin$lsociety)/sd(bos_fin$lsociety))
bos_fin$lperson_sd <- bos_fin$lperson - (mean(bos_fin$lperson)/sd(bos_fin$lperson))
bos_fin$lproperty_sd <- bos_fin$lproperty - (mean(bos_fin$lproperty)/sd(bos_fin$lproperty))

# rename racial and economic boundary variable
bos_fin$bos_white_blv <- bos_fin$p_race_white_blv
bos_fin$bos_ses_blv <- bos_fin$ses_blv

## Load Chicago data
load("chi_final.RData")

# log and +1 to variables (pop already logged)
chi_fin = chi_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime+1),
         lviolentcrime = log(violent_crime+1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_fin$larrest_sd <- chi_fin$larrests - (mean(chi_fin$larrests)/sd(chi_fin$larrests))
chi_fin$lmisdemeanors_sd <- chi_fin$lmisdemeanors - (mean(chi_fin$lmisdemeanors)/sd(chi_fin$lmisdemeanors))
chi_fin$lnonviolent_sd <- chi_fin$lnonviolent - (mean(chi_fin$lnonviolent)/sd(chi_fin$lnonviolent))
chi_fin$lsociety_sd <- chi_fin$lsociety - (mean(chi_fin$lsociety)/sd(chi_fin$lsociety))
chi_fin$lfelonies_sd <- chi_fin$lfelonies - (mean(chi_fin$lfelonies)/sd(chi_fin$lfelonies))
chi_fin$lperson_sd <- chi_fin$lperson - (mean(chi_fin$lperson)/sd(chi_fin$lperson))
chi_fin$lproperty_sd <- chi_fin$lproperty - (mean(chi_fin$lproperty)/sd(chi_fin$lproperty))
chi_fin$lviolent_sd <- chi_fin$lviolent - (mean(chi_fin$lviolent)/sd(chi_fin$lviolent))


# rename racial and economic boundary variable
chi_fin$chi_white_blv <- chi_fin$p_race_white_blv
chi_fin$chi_ses_blv <- chi_fin$ses_blv

## Load Louisville data
load("lou_final.RData")


# log and +1 to variables (pop already logged) (no misdemeanors vs felonies for louisville)
lou_fin = lou_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property+1),
         lviolentcrime = log(crime_violent+1),
         lviolent = log(violent_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
lou_fin$larrest_sd <- lou_fin$larrests - (mean(lou_fin$larrests)/sd(lou_fin$larrests))
lou_fin$lnonviolent_sd <- lou_fin$lnonviolent - (mean(lou_fin$lnonviolent)/sd(lou_fin$lnonviolent))
lou_fin$lsociety_sd <- lou_fin$lsociety - (mean(lou_fin$lsociety)/sd(lou_fin$lsociety))
lou_fin$lviolent_sd <- lou_fin$lviolent - (mean(lou_fin$lviolent)/sd(lou_fin$lviolent))
lou_fin$lperson_sd <- lou_fin$lperson - (mean(lou_fin$lperson)/sd(lou_fin$lperson))
lou_fin$lproperty_sd <- lou_fin$lproperty - (mean(lou_fin$lproperty)/sd(lou_fin$lproperty))

# rename racial and economic boundary variable
lou_fin$lou_white_blv <- lou_fin$p_race_white_blv
lou_fin$lou_ses_blv <- lou_fin$ses_blv

## Load Milwaukee data
load("mil_final.RData")


# log and +1 to variables (pop already logged)
mil_fin = mil_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_fin$larrest_sd <- mil_fin$larrests - (mean(mil_fin$larrests)/sd(mil_fin$larrests))
mil_fin$lmisdemeanors_sd <- mil_fin$lmisdemeanors - (mean(mil_fin$lmisdemeanors)/sd(mil_fin$lmisdemeanors))
mil_fin$lfelonies_sd <- mil_fin$lfelonies - (mean(mil_fin$lfelonies)/sd(mil_fin$lfelonies))
mil_fin$lnonviolent_sd <- mil_fin$lnonviolent - (mean(mil_fin$lnonviolent)/sd(mil_fin$lnonviolent))
mil_fin$lviolent_sd <- mil_fin$lviolent - (mean(mil_fin$lviolent)/sd(mil_fin$lviolent))
mil_fin$lsociety_sd <- mil_fin$lsociety - (mean(mil_fin$lsociety)/sd(mil_fin$lsociety))
mil_fin$lperson_sd <- mil_fin$lperson - (mean(mil_fin$lperson)/sd(mil_fin$lperson))
mil_fin$lproperty_sd <- mil_fin$lproperty - (mean(mil_fin$lproperty)/sd(mil_fin$lproperty))

# rename racial and economic boundary variable
mil_fin$mil_white_blv <- mil_fin$p_race_white_blv
mil_fin$mil_ses_blv <- mil_fin$ses_blv


## Load Seattle data
load("sea_final.RData")

# log and +1 to variables (pop already logged)
sea_fin = sea_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
sea_fin$larrest_sd <- sea_fin$larrests - (mean(sea_fin$larrests)/sd(sea_fin$larrests))
sea_fin$lmisdemeanors_sd <- sea_fin$lmisdemeanors - (mean(sea_fin$lmisdemeanors)/sd(sea_fin$lmisdemeanors))
sea_fin$lfelonies_sd <- sea_fin$lfelonies - (mean(sea_fin$lfelonies)/sd(sea_fin$lfelonies))
sea_fin$lnonviolent_sd <- sea_fin$lnonviolent - (mean(sea_fin$lnonviolent)/sd(sea_fin$lnonviolent))
sea_fin$lviolent_sd <- sea_fin$lviolent - (mean(sea_fin$lviolent)/sd(sea_fin$lviolent))
sea_fin$lsociety_sd <- sea_fin$lsociety - (mean(sea_fin$lsociety)/sd(sea_fin$lsociety))
sea_fin$lperson_sd <- sea_fin$lperson - (mean(sea_fin$lperson)/sd(sea_fin$lperson))
sea_fin$lproperty_sd <- sea_fin$lproperty - (mean(sea_fin$lproperty)/sd(sea_fin$lproperty))

# rename racial and economic boundary variable
sea_fin$sea_white_blv <- sea_fin$p_race_white_blv
sea_fin$sea_ses_blv <- sea_fin$ses_blv


#### reg tables associated with boundaryXcrime analysis for misdemeanor/felony arrests ####
# first create boundary dummy variable for all cities
atl_fin <- atl_fin %>% mutate(boundary_quart = quantile(atl_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

aus_fin <- aus_fin %>% mutate(boundary_quart = quantile(aus_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

bos_fin <- bos_fin %>% mutate(boundary_quart = quantile(bos_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

chi_fin <- chi_fin %>% mutate(boundary_quart = quantile(chi_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

lou_fin <- lou_fin %>% mutate(boundary_quart = quantile(lou_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

mil_fin <- mil_fin %>% mutate(boundary_quart = quantile(mil_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

sea_fin <- sea_fin %>% mutate(boundary_quart = quantile(sea_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

# run models w boundary_dummyXcrime for misdemeanors

atl.misd.crime = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus.misd.crime = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos.misd.crime = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.misd.crime = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.misd.crime = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.misd.crime = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(sea.misd.crime)


# create latex table for boundaryXcrime for misdemeanors
texreg(l = list(atl.misd.crime, aus.misd.crime, bos.misd.crime, chi.misd.crime, mil.misd.crime, sea.misd.crime),
       file = "regtable.crime2.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "lcrime" = "Log(Crime)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "boundary_quart_dummy:lcrime" = "BoundaryXCrime",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Logged Crime on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status.",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Misdemeanor Arrests" = 1:6),
       label = "table:misd_crimeXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))


# run models with boundary_dummyXcrime for felonies

atl.fel.crime = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                            phys_bound + com_density, 
                          data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus.fel.crime = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol +
                            phys_bound + com_density, 
                          data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos.fel.crime = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                            phys_bound + com_density, 
                          data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.fel.crime = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                            phys_bound + com_density, 
                          data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.fel.crime = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                            phys_bound + com_density, 
                          data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.fel.crime = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                            phys_bound + com_density, 
                          data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


# create latex table for boundaryXcrime for felonies
texreg(l = list(atl.fel.crime, aus.fel.crime, bos.fel.crime, chi.fel.crime, mil.fel.crime, sea.fel.crime),
       file = "regtable.crime1.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "lcrime" = "Log(Crime)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "boundary_quart_dummy:lcrime" = "BoundaryXCrime",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Logged Crime on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status.",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Felony Arrests" = 1:6),
       label = "table:fel_crimeXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))

lou_fin <- lou_fin %>% mutate(boundary_quart = quantile(lou_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

#### make plot for Louisville boundaryXcrime for type of arrest

# Arrests against persons
# expected values
mdata <- lou_fin %>% dplyr::select(lperson_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()

lou1 = lm_robust(lperson_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                   phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,11.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(lou1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,11.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(lou1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

lou_preds_person <- preds %>% mutate(type = "Against Persons")


# Arrests against property
# expected values
mdata <- lou_fin %>% dplyr::select(lproperty_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()

lou1 = lm_robust(lproperty_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                   phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,11.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(lou1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,11.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(lou1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

lou_preds_property <- preds %>% mutate(type = "Against Property")

# Arrests against society
# expected values
mdata <- lou_fin %>% dplyr::select(lsociety_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()

lou1 = lm_robust(lsociety_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,11.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(lou1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,11.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(lou1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

lou_preds_society <- preds %>% mutate(type = "Against Society")

## put together all types and graph)
## expected values graph ##
lou_preds_person <- lou_preds_person %>% dplyr::select(lcrime,fitted.fit,
                                                       fitted.lwr,fitted.upr,boundary,type)
lou_preds_property <- lou_preds_property %>% dplyr::select(lcrime,fitted.fit,
                                                           fitted.lwr,fitted.upr,boundary,type)
lou_preds_society <- lou_preds_society %>% dplyr::select(lcrime,fitted.fit,
                                                         fitted.lwr,fitted.upr,boundary,type)



d <- rbind(lou_preds_person, lou_preds_property, lou_preds_society)

head(d)

p <- ggplot(data=d, aes(x=lcrime, y=fitted.fit,group=boundary)) +
  facet_wrap(~ type, nrow = 3, scales = "free") +
  #geom_point(aes(shape = Boundary), color = "grey",alpha=.5)+
  geom_line(aes(y=fitted.fit, linetype=boundary, color = boundary))+
  geom_ribbon(aes(ymin=fitted.lwr, ymax=fitted.upr,linetype=boundary, color=boundary),alpha=0.1)+
  labs(y="Arrests (Logged)",x="Crime (Logged)")+
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA))+
  theme(strip.background = element_rect(fill="white"))


ggsave(plot = p, filename = "louisville_raceboundXcrime_outsample.png", 
       height = 8, width = 8)

#### reg tables associated with boundaryXcrime analysis for against person/property/society arrests ####

# against person


lou.pers.crime = lm_robust(lperson_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

lou.prop.crime = lm_robust(lproperty_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                             phys_bound + com_density, 
                           data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

lou.soc.crime = lm_robust(lsociety_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + 
                            phys_bound + com_density, 
                          data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


# create latex table for boundaryXcrime for against persons
texreg(l = list(lou.pers.crime, lou.prop.crime, lou.soc.crime),
       file = "regtable.crimelou.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "lcrime" = "Log(Crime)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "boundary_quart_dummy:lcrime" = "BoundaryXCrime",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Lousiville: Influence of Logged Crime on Logged Arrests by Type (Standardized), Conditional on Racial Boundary Status.",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("Against Persons","Against Property","Against Society"),
       label = "table:lou_crimeXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))



